PR-S003-S015_data-analysis
Jenna Hanlon
2024-01-11
Load Libraries
library(qiime2R)
library(ggplot2)
library(vegan) ## Warning: package 'vegan' was built under R version 4.1.3## Loading required package: permute## Warning: package 'permute' was built under R version 4.1.3## Loading required package: lattice## Warning: package 'lattice' was built under R version 4.1.3## This is vegan 2.6-4library(plyr) ## Warning: package 'plyr' was built under R version 4.1.3library(dplyr) ## Warning: package 'dplyr' was built under R version 4.1.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(scales)## Warning: package 'scales' was built under R version 4.1.3library(grid)
library(reshape2)## Warning: package 'reshape2' was built under R version 4.1.3library(phyloseq)
library(picante) ## Warning: package 'picante' was built under R version 4.1.3## Loading required package: ape## Warning: package 'ape' was built under R version 4.1.3##
## Attaching package: 'ape'## The following object is masked from 'package:dplyr':
##
## where## Loading required package: nlme## Warning: package 'nlme' was built under R version 4.1.3##
## Attaching package: 'nlme'## The following object is masked from 'package:dplyr':
##
## collapselibrary(tidyr) ## Warning: package 'tidyr' was built under R version 4.1.3##
## Attaching package: 'tidyr'## The following object is masked from 'package:reshape2':
##
## smithslibrary(viridis) ## Loading required package: viridisLite##
## Attaching package: 'viridis'## The following object is masked from 'package:scales':
##
## viridis_pallibrary(qiime2R)
library(DESeq2) ## Loading required package: S4Vectors## Warning: package 'S4Vectors' was built under R version 4.1.3## Loading required package: stats4## Loading required package: BiocGenerics##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:dplyr':
##
## combine, intersect, setdiff, union## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, basename, cbind, colnames,
## dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
## grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
## rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
## union, unique, unsplit, which.max, which.min##
## Attaching package: 'S4Vectors'## The following object is masked from 'package:tidyr':
##
## expand## The following objects are masked from 'package:dplyr':
##
## first, rename## The following object is masked from 'package:plyr':
##
## rename## The following objects are masked from 'package:base':
##
## expand.grid, I, unname## Loading required package: IRanges##
## Attaching package: 'IRanges'## The following object is masked from 'package:nlme':
##
## collapse## The following object is masked from 'package:phyloseq':
##
## distance## The following objects are masked from 'package:dplyr':
##
## collapse, desc, slice## The following object is masked from 'package:plyr':
##
## desc## The following object is masked from 'package:grDevices':
##
## windows## Loading required package: GenomicRanges## Loading required package: GenomeInfoDb## Loading required package: SummarizedExperiment## Loading required package: MatrixGenerics## Loading required package: matrixStats## Warning: package 'matrixStats' was built under R version 4.1.3##
## Attaching package: 'matrixStats'## The following object is masked from 'package:dplyr':
##
## count## The following object is masked from 'package:plyr':
##
## count##
## Attaching package: 'MatrixGenerics'## The following objects are masked from 'package:matrixStats':
##
## colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
## colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
## colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
## colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
## colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
## colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
## colWeightedMeans, colWeightedMedians, colWeightedSds,
## colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
## rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
## rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
## rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
## rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
## rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
## rowWeightedSds, rowWeightedVars## Loading required package: Biobase## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.##
## Attaching package: 'Biobase'## The following object is masked from 'package:MatrixGenerics':
##
## rowMedians## The following objects are masked from 'package:matrixStats':
##
## anyMissing, rowMedians## The following object is masked from 'package:phyloseq':
##
## sampleNameslibrary(patchwork)
library(RColorBrewer)## Warning: package 'RColorBrewer' was built under R version 4.1.3library(microViz) ## microViz version 0.10.10 - Copyright (C) 2023 David Barnett
## ! Website: https://david-barnett.github.io/microViz
## v Useful? For citation details, run: `citation("microViz")`
## x Silence? `suppressPackageStartupMessages(library(microViz))`library(speedyseq) ##
## Attaching package: 'speedyseq'## The following objects are masked from 'package:phyloseq':
##
## filter_taxa, plot_bar, plot_heatmap, plot_tree, psmelt, tax_glom,
## tip_glom, transform_sample_countslibrary(ComplexHeatmap) ## ========================================
## ComplexHeatmap version 2.10.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
##
## If you use it in published research, please cite:
## Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
## genomic data. Bioinformatics 2016.
##
## The new InteractiveComplexHeatmap package can directly export static
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================library(ggVennDiagram) ## Warning: package 'ggVennDiagram' was built under R version 4.1.3library(SuperExactTest) ##
## Attaching package: 'SuperExactTest'## The following objects are masked from 'package:GenomicRanges':
##
## intersect, union## The following object is masked from 'package:GenomeInfoDb':
##
## intersect## The following objects are masked from 'package:IRanges':
##
## intersect, union## The following objects are masked from 'package:S4Vectors':
##
## intersect, union## The following objects are masked from 'package:BiocGenerics':
##
## intersect, union## The following objects are masked from 'package:dplyr':
##
## intersect, union## The following objects are masked from 'package:base':
##
## intersect, unionlibrary(nVennR)
library(vegan)
library(plyr)
library(dplyr)
library(scales)
library(grid)
library(reshape2)
library(phyloseq)
library(picante)
library(tidyr)
library(viridis)
library(qiime2R)
library(DESeq2)
library(patchwork)
library(RColorBrewer)
library(speedyseq)
library(RColorBrewer)
library(microViz)
library(ComplexHeatmap)
library(plotly) ##
## Attaching package: 'plotly'## The following object is masked from 'package:ComplexHeatmap':
##
## add_heatmap## The following object is masked from 'package:microViz':
##
## add_paths## The following object is masked from 'package:IRanges':
##
## slice## The following object is masked from 'package:S4Vectors':
##
## rename## The following objects are masked from 'package:plyr':
##
## arrange, mutate, rename, summarise## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutlibrary(seecolor)
library(htmlwidgets)## Warning: package 'htmlwidgets' was built under R version 4.1.3library(htmltools)## Warning: package 'htmltools' was built under R version 4.1.3set.seed(17384)Loading Data
PRphyseq <- qza_to_phyloseq(features="S003-S015-table.qza", tree="S003-S015-rooted-tree.qza",taxonomy="S003-S015-taxonomy.qza", metadata = "S003-S015-Field-Measurements-2.txt")rank_names(PRphyseq)## [1] "Kingdom" "Phylum" "Class" "Order" "Family" "Genus" "Species"sample_variables(PRphyseq)## [1] "Project" "Sub_Location"
## [3] "Location" "Season"
## [5] "Month" "Year"
## [7] "Sample_Type" "Sample_ID"
## [9] "Day" "Station_ID"
## [11] "Station_Name" "Collected_by"
## [13] "Time" "Samples_Collected_by"
## [15] "Temperature_C" "mm_Hg"
## [17] "DO_mg_L" "DO_LDO"
## [19] "Specific_Conductance_uS_cm" "Salinity_ppt"
## [21] "pH" "pH_mv"
## [23] "IVCH_ug_L" "Turbidity_.Fluorometer"
## [25] "Comments"Removing chloroplast sequences and any contaminant sequences
## Removing chloroplast sequences and any contaminant sequences
PRphyseq <- subset_taxa(PRphyseq, Kingdom != "d__Eukaryota")
PRphyseq <- subset_taxa(PRphyseq, Kingdom != "d__Archaea")
PRphyseq <- subset_taxa(PRphyseq, Order != "o__Chloroplast")
PRphyseq <- subset_taxa(PRphyseq, Order != "Chloroplast")
PRphyseq <- subset_taxa(PRphyseq, Family != "f__Mitochondria")
PRphyseq <- subset_taxa(PRphyseq, Family != "Mitochondria")
PRphyseq <- subset_taxa(PRphyseq, Family != "Chloroplast")
##Check that contaminant sequences are removed (easiest to save as data frame and search) taxtabl <- as.data.frame(tax_table(S003physeq)) #At this point, blanks and positive controls should also be removed #S003physeq = subset_samples(S003physeq, Sample != "mock")View((tax_table(PRphyseq)))
##This will allow you to open your data frame. You can search in the top right corner of the dataframe for the sequences you wanted to remove in above code to ensure they were removed. Creating a subset of samples
PRphyseq <- PRphyseq %>% subset_samples(Project %in% c("S003","S015")) ##This will create a subset (based on project column) of the samples we would like to use in the analysis. If you would like to keep all samples in the metadata file, you can ignore this step.
PR_Normal <- PRphyseq %>% subset_samples(Sample_Type %in% c("N"))Examine the number of reads
readsumsdf = data.frame(nreads = sort(taxa_sums(PR_Normal), TRUE), sorted = 1:ntaxa(PR_Normal), type = "ASVs")
readsumsdf = rbind(readsumsdf, data.frame(nreads = sort(sample_sums(PR_Normal), TRUE), sorted = 1:nsamples(PR_Normal), type = "Samples"))
title = "Total number of reads"
p = ggplot(readsumsdf, aes(x = sorted, y = nreads)) + geom_bar(stat = "identity")
p + ggtitle(title) + scale_y_log10() + facet_wrap(~type, 1, scales = "free")## Warning: Transformation introduced infinite values in continuous y-axis## Warning: Removed 1908 rows containing missing values (`geom_bar()`).
Plot a rarefaction curve using vegan:
taxa_are_rows(PR_Normal)## [1] TRUEmat <- t(otu_table(PR_Normal))
class(mat) <- "matrix"## Warning in class(mat) <- "matrix": Setting class(x) to "matrix" sets attribute
## to NULL; result will no longer be an S4 objectclass(mat)## [1] "matrix" "array"mat <- as(t(otu_table(PR_Normal)), "matrix")
class(mat)## [1] "matrix" "array"raremax <- min(rowSums(mat))
system.time(rarecurve(mat, step = 100, sample = raremax, col = "blue", label = FALSE))
## user system elapsed
## 4.28 0.17 6.32Creating table for number of reads per sample
#Sum reads and sort by samples with least to greatest number of reads
samples_reads <- sort(sample_sums(PR_Normal))
#Create table
knitr::kable(samples_reads) %>%
kableExtra::kable_styling("striped", latex_options="scale_down") %>%
kableExtra::scroll_box(width = "100%")| x | |
|---|---|
| PR109C | 14 |
| pr242a | 679 |
| PR325 | 699 |
| PR164C | 700 |
| PR151C | 855 |
| PR566A | 1040 |
| PR321 | 1088 |
| PR147C | 1128 |
| PR564A | 1219 |
| lsj3-1 | 1251 |
| PR323 | 1323 |
| PR559A | 1335 |
| PR145C | 1351 |
| PR422 | 1363 |
| PR425 | 1372 |
| PR162C | 1450 |
| PR141C | 1461 |
| PR421 | 1479 |
| PR317 | 1481 |
| PR399 | 1579 |
| PR561A | 1624 |
| PR311 | 1634 |
| PR420 | 1700 |
| PR139C | 1730 |
| PR333 | 1785 |
| PR558A | 1793 |
| PR554A | 1972 |
| PR419 | 2037 |
| PR417 | 2075 |
| PR160C | 2191 |
| PR107C | 2201 |
| PR556A | 2311 |
| PR334 | 2388 |
| PR557A | 2532 |
| PR392 | 2952 |
| PR395 | 3001 |
| PR149C | 3054 |
| PR394 | 3055 |
| PR332 | 3152 |
| PR393 | 3167 |
| PR555A | 3238 |
| PR416 | 3539 |
| PR562A | 3614 |
| PR415 | 3728 |
| PR328 | 3792 |
| PR400 | 3815 |
| PR316 | 4073 |
| PR565A | 4122 |
| PR331 | 4337 |
| PR402 | 4429 |
| PR560A | 4453 |
| PR401 | 4461 |
| PR563A | 4999 |
| PR315 | 5085 |
| 239 | 5356 |
| 18 | 5464 |
| pr248a | 5585 |
| PR318 | 5696 |
| PR102C | 6357 |
| PR322 | 6779 |
| pr250a | 7781 |
| pr262a | 8718 |
| pr271a | 9101 |
| 101 | 9859 |
| pr267a | 10277 |
| 220 | 10514 |
| pr247a | 10711 |
| 44 | 10755 |
| 104 | 10896 |
| cs1-1 | 11179 |
| cmpw-3 | 11765 |
| pr295a | 12285 |
| lsj3-2 | 12316 |
| cmpe-3 | 12325 |
| pr249a | 12448 |
| cmp2-2 | 12550 |
| 142 | 12638 |
| 198 | 12984 |
| lp1-2 | 13656 |
| 42 | 13896 |
| cmp2-1 | 13943 |
| cmpe-2 | 14167 |
| cmp2-3 | 14339 |
| lp1-3 | 14775 |
| 5 | 14776 |
| 218 | 15212 |
| 181 | 15472 |
| cmpe-1 | 15558 |
| cmpw-1 | 15630 |
| 183 | 16414 |
| 63 | 16771 |
| lsj3-3 | 16802 |
| cs1-2 | 19707 |
| 95 | 21854 |
| lp1-1 | 23460 |
| 106 | 25062 |
| 90 | 25753 |
| cs1-3 | 29344 |
| cmpw-2 | 51074 |
| pr293a | 73188 |
Remove samples with <1,000 reads
#Remove samples with less than 1,000 reads
PR_Normal = subset_samples(PR_Normal,
Sample_ID != "PR105C")
PR_Normal = subset_samples(PR_Normal,
Sample_ID != "PR109C")
PR_Normal = subset_samples(PR_Normal,
Sample_ID != "pr242a")
PR_Normal = subset_samples(PR_Normal,
Sample_ID != "PR325")
PR_Normal = subset_samples(PR_Normal,
Sample_ID != "PR164C")
PR_Normal = subset_samples(PR_Normal,
Sample_ID != "PR151C")Examining abundance of individual ASVs
#For each ASV, find the number of samples in which each ASV is 0, then divide by the total number of samples
test_filter_method <- as.data.frame(rowSums(otu_table(PR_Normal) == 0)/ncol(otu_table(PR_Normal)))
#Change the name of the column in the test_filter_method4 dataframe; this column contains the proportion of samples in which each ASV is 0
colnames(test_filter_method) <- c('samplew0')
#Create a bar plot
ggplot(data = test_filter_method) + geom_bar(mapping = aes(x= samplew0)) +labs(x = "Proportion of samples in which ASV = 0", y = "# of ASV's") +theme(text = element_text(size = 18), axis.title = element_text(size = 15),panel.spacing = unit(1, "lines"),panel.border = element_rect(colour = "black", fill = NA, size = 0.5), panel.background = element_blank()) ## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## i Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Removing singletons and doubletons
##How many singletons and doubletons are present
#Create a dataframe that includes only the total number of times that each ASV occurs
readsumsdf_ASVs <- readsumsdf %>%
filter(type == "ASVs")
#How many singletons are present?
length(which(readsumsdf_ASVs$nreads <= 0))## [1] 1908length(which(readsumsdf_ASVs$nreads == 1))## [1] 7#How many doubletons are present?
length(which(readsumsdf_ASVs$nreads == 2))## [1] 158#Filter out low abundance ASVs
PR_Normal = prune_taxa(taxa_sums(PR_Normal) > 2, PR_Normal)Examine the abundance of individual ASVs again:
#For each ASV, find the number of samples in which each ASV is 0, then divide by the total number of samples
test_filter_method_filtered <- as.data.frame(rowSums(otu_table(PR_Normal) == 0)/ncol(otu_table(PR_Normal)))
#Change the name of the column in the test_filter_method4 dataframe; this column contains the proportion of samples in which each ASV is 0
colnames(test_filter_method_filtered) <- c('samplew0')
#Create a bar plot
ggplot(data = test_filter_method_filtered) + geom_bar(mapping = aes(x= samplew0)) +
labs(x = "Proportion of samples in which ASV = 0", y = "# of ASV's") +
theme(text = element_text(size = 18), axis.title = element_text(size = 15),
panel.spacing = unit(1, "lines"),
panel.border = element_rect(colour = "black", fill = NA, size = 0.5),
panel.background = element_blank()) 
Label variables
sample_data(PR_Normal)$Sub_Location <- factor(sample_data(PR_Normal)$Sub_Location, levels = c("BSJ1","BSJ1-ENR","BSJ2","BSJ2-ENR","BSJ3","BSJ3-ENR","BSJ4-ENR","BSJ5-ENR","BSJ6-ENR","BSJ7-ENR","CMP1","CMP1-ENR","CMP2","CMP2-ENR","CMPE","CMPW","CS1","CS1-ENR","CS2","CS2-ENR","CSA1-ENR","LLC1","LLC1-ENR","LP1","LP1-ENR","LP1-ENR","LP1L-ENR","LP2-ENR","LP3-ENR","LP4-ENR","LP5-ENR","LSJ1","LSJ1-ENR","LSJ2","LSJ2-ENR","LSJ3","LT1","LT1-ENR","LT2","LT2-ENR","LT3-ENR","QB1-ENR","QSA1-ENR","RPN1-ENR"),
labels = c("BSJ1","BSJ1-ENR","BSJ2","BSJ2-ENR","BSJ3","BSJ3-ENR","BSJ4-ENR","BSJ5-ENR","BSJ6-ENR","BSJ7-ENR","CMP1","CMP1-ENR","CMP2","CMP2-ENR","CMPE","CMPW","CS1","CS1-ENR","CS2","CS2-ENR","CSA1-ENR","LLC1","LLC1-ENR","LP1","LP1-ENR","LP1-ENR","LP1L-ENR","LP2-ENR","LP3-ENR","LP4-ENR","LP5-ENR","LSJ1","LSJ1-ENR","LSJ2","LSJ2-ENR","LSJ3","LT1","LT1-ENR","LT2","LT2-ENR","LT3-ENR","QB1-ENR","QSA1-ENR","RPN1-ENR"))
sample_data(PR_Normal)$Location <- factor(sample_data(PR_Normal)$Location, levels = c(
"BSJ","BSJ-ENR","CMP","CMP-ENR","CMPE","CMPW", "CS","CS-ENR", "CSA","CSA-ENR","LLC","LLC-ENR","LP","LP-ENR","LSJ","LSJ-ENR","LT","LT-ENR","QB","QB-ENR","QSA","QSA-ENR","RPN","RPN-ENR"),
labels = c("BSJ","BSJ-ENR","CMP","CMP-ENR","CMPE","CMPW", "CS","CS-ENR", "CSA","CSA-ENR","LLC","LLC-ENR","LP","LP-ENR","LSJ","LSJ-ENR","LT","LT-ENR","QB","QB-ENR","QSA","QSA-ENR","RPN","RPN-ENR"))
sample_data(PR_Normal)$Season <- factor(sample_data(PR_Normal)$Season, levels = c("W", "D"),
labels = c("W", "D"))
sample_data(PR_Normal)$Project <- factor(sample_data(PR_Normal)$Project,
levels = c("S003","S015"), labels = c("S003","S015"))
sample_data(PR_Normal)$Month <- factor(sample_data(PR_Normal)$Month,
levels = c("April","August","December","January","July","June","March","May","November","October","September"), labels = c("April","August","December","January","July","June","March","May","November","October","September"))
sample_data(PR_Normal)$Year <- factor(sample_data(PR_Normal)$Year,
levels = c("2021","2022","2023"), labels = c("2021","2022","2023"))
sample_data(PR_Normal)$Sample_Type <- factor(sample_data(PR_Normal)$Sample_Type,
levels = c("N","E"), labels = c("N","E"))
sample_data(PR_Normal)$Sample_ID <- factor(sample_data(PR_Normal)$Sample_ID,
levels = c("PR419","PR417","PR422","PR421","PR420","PR416","PR415","90","44","42","63","5","18","pr293a", "pr271a","pr262a","pr285a","pr295a","pr267a","pr290a","PR425","PR325","PR323","PR321","PR322","PR332","PR333","PR331","PR315","PR318","PR328","PR317","PR316","PR334", "PR311","PR563A","PR562A","PR560A","PR561A","PR565A","cmp2-3","cmpe-3","cmpw-3","cs1-3","PR556A","PR559A","PR564A","lp1-3","PR558A","PR557A","PR566A","lsj3-3","PR555A","PR554A","cmp2-2","cmpe-2","cmpw-2","cs1-2","lp1-2","lsj3-2","PR395","PR394","PR392", "PR393","PR400","PR401","PR399","PR402","PR464A","PR465A","PR466A","PR467A","PR468A","PR469A","PR470A","PR473A","cmp2-1","PR472A","cmpe-1","cmpw-1","cs1-1","PR476A","PR477A","PR474A","PR478A","lp1-1","PR488A","PR489A","PR490A","PR491A","PR492A","PR493A","PR479A","PR480A","lsj3-1","PR481A","PR482A","PR483A","PR484A","PR485A","PR471A","Mock","Undeter","B53","B84","B59","B43","B36","B44","B107","B90","B126","B141","B128","B152","B134","B132","B140","B89","B100","B102","B60","B76","B56","blank","blank2","Mock-com","B130","B136","B138","B33","B25","B41","B81","B68","B57","B87","B91","B117","B143","B139","B154","B69","B52","B75","B99","B118","B124","B26","B34","B42","pr249a","pr248a","pr250a","pr242a","pr247a","183","218","181","198","239","220","PR160C","PR162C","PR164C","PR102C","101","PR109C","108","PR107C","106","PR143C","142","PR145C","PR151C","PR96C","95","PR147C","PR149C","PR105C","104","PR139C","PR141C"
),
labels = c("PR419","PR417","PR422","PR421","PR420","PR416","PR415","90","44","42","63","5","18","pr293a", "pr271a","pr262a","pr285a","pr295a","pr267a","pr290a","PR425","PR325","PR323","PR321","PR322","PR332","PR333","PR331","PR315","PR318","PR328","PR317","PR316","PR334", "PR311","PR563A","PR562A","PR560A","PR561A","PR565A","cmp2-3","cmpe-3","cmpw-3","cs1-3","PR556A","PR559A","PR564A","lp1-3","PR558A","PR557A","PR566A","lsj3-3","PR555A","PR554A","cmp2-2","cmpe-2","cmpw-2","cs1-2","lp1-2","lsj3-2","PR395","PR394","PR392", "PR393","PR400","PR401","PR399","PR402","PR464A","PR465A","PR466A","PR467A","PR468A","PR469A","PR470A","PR473A","cmp2-1","PR472A","cmpe-1","cmpw-1","cs1-1","PR476A","PR477A","PR474A","PR478A","lp1-1","PR488A","PR489A","PR490A","PR491A","PR492A","PR493A","PR479A","PR480A","lsj3-1","PR481A","PR482A","PR483A","PR484A","PR485A","PR471A","Mock","Undeter","B53","B84","B59","B43","B36","B44","B107","B90","B126","B141","B128","B152","B134","B132","B140","B89","B100","B102","B60","B76","B56","blank","blank2","Mock-com","B130","B136","B138","B33","B25","B41","B81","B68","B57","B87","B91","B117","B143","B139","B154","B69","B52","B75","B99","B118","B124","B26","B34","B42","pr249a","pr248a","pr250a","pr242a","pr247a","183","218","181","198","239","220","PR160C","PR162C","PR164C","PR102C","101","PR109C","108","PR107C","106","PR143C","142","PR145C","PR151C","PR96C","95","PR147C","PR149C","PR105C","104","PR139C","PR141C"))Tax Fix
#View the tax_table
knitr::kable(head(tax_table(PR_Normal))) %>%
kableExtra::kable_styling("striped") %>%
kableExtra::scroll_box(width = "100%") | Kingdom | Phylum | Class | Order | Family | Genus | Species | |
|---|---|---|---|---|---|---|---|
| 1779c60c03d948a73687f440be973f68 | d__Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Flammeovirgaceae | Algivirga | Algivirga_pacifica |
| 25feebceed837c1dfa810cd8c0f26b41 | d__Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Cyclobacteriaceae | NA | NA |
| 0872b535e6e15b9cd031dae6a67d6257 | d__Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Amoebophilaceae | uncultured | uncultured_bacterium |
| 6590e42d5aeb18c338a836bc26028592 | d__Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | uncultured | uncultured | uncultured_Bacteroidetes |
| 0f22f07d74171ba4bbe6b3977217cbdc | d__Bacteria | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidetes_BD2-2 | Bacteroidetes_BD2-2 | uncultured_organism |
| d817bb1b166e13e13db6db8d56b43983 | d__Bacteria | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidetes_BD2-2 | Bacteroidetes_BD2-2 | NA |
#Now fix the labels
PR_Normal <- tax_fix(PR_Normal)
#View the relabeled tax_table
knitr::kable(head(tax_table(PR_Normal))) %>%
kableExtra::kable_styling("striped") %>%
kableExtra::scroll_box(width = "100%")| Kingdom | Phylum | Class | Order | Family | Genus | Species | |
|---|---|---|---|---|---|---|---|
| 1779c60c03d948a73687f440be973f68 | d__Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Flammeovirgaceae | Algivirga | Algivirga_pacifica |
| 25feebceed837c1dfa810cd8c0f26b41 | d__Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Cyclobacteriaceae | Cyclobacteriaceae Family | Cyclobacteriaceae Family |
| 0872b535e6e15b9cd031dae6a67d6257 | d__Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Amoebophilaceae | uncultured | uncultured_bacterium |
| 6590e42d5aeb18c338a836bc26028592 | d__Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | uncultured | uncultured | uncultured_Bacteroidetes |
| 0f22f07d74171ba4bbe6b3977217cbdc | d__Bacteria | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidetes_BD2-2 | Bacteroidetes_BD2-2 | uncultured_organism |
| d817bb1b166e13e13db6db8d56b43983 | d__Bacteria | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidetes_BD2-2 | Bacteroidetes_BD2-2 | Bacteroidetes_BD2-2 Genus |
Bacterial community composition bar plots
library(seecolor)
bcc_hex <- print_color(c("black", "#440154FF", "#450659FF","#460B5EFF","#472D7AFF","#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF", "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF","#2E6E8EFF", "#2D718EFF", "#2B748EFF","#2A778EFF", "#297B8EFF","#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF","#1FA287FF", "#21A585FF", "#23A983FF","#25AC82FF", "#29AF7FFF", "#2DB27DFF","#32B67AFF", "#37B878FF", "#3CBC74FF","#57C766FF", "#5EC962FF","#A2DA37FF","#DAE319FF", "#E4E419FF","#ECE51BFF", "#F5E61FFF","darkorchid1", "darkorchid2", "darkorchid3","#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF","#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#FF6699","#F68D45FF", "#FCA537FF","#F6E726FF", "#F4ED27FF", "#00489C", "#CCCCCC", "#999999", "#A1C299","#300018"), type = "r")##
## ------ c black #440154FF #450659FF #460B5EFF #472D7AFF #3B518BFF #3A548CFF #38598CFF #365C8DFF #34608DFF #33638DFF #31678EFF #306A8EFF #2E6E8EFF #2D718EFF #2B748EFF #2A778EFF #297B8EFF #1F958BFF #1F988BFF #1E9C89FF #1F9F88FF #1FA287FF #21A585FF #23A983FF #25AC82FF #29AF7FFF #2DB27DFF #32B67AFF #37B878FF #3CBC74FF #57C766FF #5EC962FF #A2DA37FF #DAE319FF #E4E419FF #ECE51BFF #F5E61FFF darkorchid1 darkorchid2 darkorchid3 #A21C9AFF #A62098FF #AB2394FF #AE2892FF #B22B8FFF #B6308BFF #BA3388FF #BE3885FF #C13B82FF #C53F7EFF #C8437BFF #CC4678FF #CE4B75FF #D14E72FF #D5536FFF #D7566CFF #DA5B69FF #DD5E66FF #E06363FF #FF6699 #F68D45FF #FCA537FF #F6E726FF #F4ED27FF #00489C #CCCCCC #999999 #A1C299 #300018 ------
## black
## #440154FF
## #450659FF
## #460B5EFF
## #472D7AFF
## #3B518BFF
## #3A548CFF
## #38598CFF
## #365C8DFF
## #34608DFF
## #33638DFF
## #31678EFF
## #306A8EFF
## #2E6E8EFF
## #2D718EFF
## #2B748EFF
## #2A778EFF
## #297B8EFF
## #1F958BFF
## #1F988BFF
## #1E9C89FF
## #1F9F88FF
## #1FA287FF
## #21A585FF
## #23A983FF
## #25AC82FF
## #29AF7FFF
## #2DB27DFF
## #32B67AFF
## #37B878FF
## #3CBC74FF
## #57C766FF
## #5EC962FF
## #A2DA37FF
## #DAE319FF
## #E4E419FF
## #ECE51BFF
## #F5E61FFF
## darkorchid1
## darkorchid2
## darkorchid3
## #A21C9AFF
## #A62098FF
## #AB2394FF
## #AE2892FF
## #B22B8FFF
## #B6308BFF
## #BA3388FF
## #BE3885FF
## #C13B82FF
## #C53F7EFF
## #C8437BFF
## #CC4678FF
## #CE4B75FF
## #D14E72FF
## #D5536FFF
## #D7566CFF
## #DA5B69FF
## #DD5E66FF
## #E06363FF
## #FF6699
## #F68D45FF
## #FCA537FF
## #F6E726FF
## #F4ED27FF
## #00489C
## #CCCCCC
## #999999
## #A1C299
## #300018Creating Family Barplot
##Subsetting non-enriched samples to view composition in bar plot
PR_N <- PR_Normal %>% subset_samples(Sample_Type %in% c("N"))PRSeqR_family <- PR_N %>%
tax_glom(taxrank = "Family") %>%
transform_sample_counts(function(x) {x/sum(x)}) %>%
psmelt() %>%
group_by(Location,
Kingdom, Phylum, Class, Order, Family) %>%
filter(Abundance > 0.01) %>%
arrange(Class)
PRSeqR_family$Class_family <- paste(PRSeqR_family$Class, PRSeqR_family$Family, sep="_")
PRSeqR_family_bar <- ggplot(PRSeqR_family, aes(x = Location, y = Abundance, fill = Class_family)) +
geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") + scale_fill_manual(values = c("black", "#440154FF", "#450659FF","#460B5EFF", "#472D7AFF", "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF", "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF", "#2E6E8EFF", "#2D718EFF", "#2B748EFF", "#2A778EFF", "#297B8EFF", "#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF", "#1FA287FF", "#21A585FF", "#23A983FF", "#25AC82FF", "#29AF7FFF", "#2DB27DFF", "#32B67AFF", "#37B878FF", "#3CBC74FF", "#57C766FF", "#5EC962FF", "#A2DA37FF", "#DAE319FF", "#E4E419FF", "#ECE51BFF", "#F5E61FFF", "darkorchid1", "darkorchid2", "darkorchid3", "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF", "#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#00CCCC", "#006600", "#FF6699", "#F68D45FF", "#FCA537FF", "#F6E726FF", "#F4ED27FF", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#191970","#353935","#F0FFFF","#89CFF0","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58")) +
guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) +
ylab("Relative Abundance (Family > 1%) \n") + xlab("Location") +
theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"), panel.grid.minor = element_blank(),
panel.grid.major = element_blank(), axis.text.x =element_text(angle = 90, vjust = 0.5, hjust=1, size=8, colour="black"),
axis.text.y = element_text(size=8, colour="black"), plot.title = element_text(hjust = 0.5),
axis.ticks.x = element_line(colour="#000000", size=0.1), axis.ticks.y = element_line(colour="#000000", size=0.1),
legend.text = element_text(size = 7))## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## i Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## i Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.print(PRSeqR_family_bar)
p <- ggplotly(PRSeqR_family_bar)htmlwidgets::saveWidget(p, "PRSeqR_family_bar.html")
htmltools::tags$iframe(
src=file.path(getwd(), "PRSeqR_family_bar.html"),
width="100%",
height="600",
scrolling="no",
seamless="seamless",
frameBorder="0")Creating Species Barplot
PRSeqR_Species <- PR_N %>%
tax_glom(taxrank = "Species") %>%
transform_sample_counts(function(x) {x/sum(x)}) %>%
psmelt() %>%
group_by(Location,
Kingdom, Phylum, Class, Order, Family,Genus,Species) %>%
filter(Abundance > 0.05) %>%
arrange(Species)
PRSeqR_Species_bar <- ggplot(PRSeqR_Species, aes(x = Location, y = Abundance, fill = Species)) +
geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") + scale_fill_manual(values = c("black", "#440154FF", "#450659FF","#460B5EFF", "#472D7AFF", "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF", "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF", "#2E6E8EFF", "#2D718EFF", "#2B748EFF", "#2A778EFF", "#297B8EFF", "#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF", "#1FA287FF", "#21A585FF", "#23A983FF", "#25AC82FF", "#29AF7FFF", "#2DB27DFF", "#32B67AFF", "#37B878FF", "#3CBC74FF", "#57C766FF", "#5EC962FF", "#A2DA37FF", "#DAE319FF", "#E4E419FF", "#ECE51BFF", "#F5E61FFF", "darkorchid1", "darkorchid2", "darkorchid3", "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF", "#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#00CCCC", "#006600", "#FF6699", "#F68D45FF", "#FCA537FF", "#F6E726FF", "#F4ED27FF", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#191970","#353935","#F0FFFF","#89CFF0","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58")) +
guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) +
ylab("Relative Abundance (Species > 5%) \n") + xlab("Sample ID") +
theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"), panel.grid.minor = element_blank(),
panel.grid.major = element_blank(), axis.text.x =element_text(angle = 90, vjust = 0.5, hjust=1, size=8, colour="black"),
axis.text.y = element_text(size=8, colour="black"), plot.title = element_text(hjust = 0.5),
axis.ticks.x = element_line(colour="#000000", size=0.1), axis.ticks.y = element_line(colour="#000000", size=0.1),
legend.text = element_text(size = 7))p2 <- ggplotly(PRSeqR_Species_bar)print(PRSeqR_Species_bar)
htmlwidgets::saveWidget(p, "PRSeqR_Species_bar.html")
htmltools::tags$iframe(
src=file.path(getwd(), "PRSeqR_Species_bar.html"),
width="100%",
height="600",
scrolling="no",
seamless="seamless",
frameBorder="0")Creating Season vs Location Barplot with non-enriched samples
PRSeqR_family_Season <- PR_N %>%
tax_glom(taxrank = "Family") %>%
transform_sample_counts(function(x) {x/sum(x)}) %>%
psmelt() %>%
group_by(Season,
Kingdom, Phylum, Class, Order, Family)%>%
dplyr::summarize(Mean =
mean(Abundance, na.rm=TRUE)) %>%
filter(Mean > 0.01) %>%
arrange(Class) ## `summarise()` has grouped output by 'Season', 'Kingdom', 'Phylum', 'Class',
## 'Order'. You can override using the `.groups` argument.## `summarise()` has grouped output by 'Season', 'Kingdom', 'Phylum', 'Class',
## 'Order'. You can override using the `.groups` argument.## Creaitng the plot
PRSeqR_family_Season$Class_family <- paste(PRSeqR_family_Season$Class, PRSeqR_family_Season$Family, sep="_")
PRSeqR_Season_family_bar <- ggplot(PRSeqR_family_Season, aes(x = Season, y = Mean, fill = Class_family)) +
geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") + scale_fill_manual(values = c("black", "#440154FF", "#450659FF","#460B5EFF","#472D7AFF", "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF","#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF","#2E6E8EFF", "#2D718EFF", "#2B748EFF","#2A778EFF", "#297B8EFF","#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF", "#1FA287FF", "#21A585FF", "#23A983FF","#25AC82FF", "#29AF7FFF", "#2DB27DFF","#32B67AFF", "#37B878FF", "#3CBC74FF","#57C766FF", "#5EC962FF","#A2DA37FF", "#DAE319FF", "#E4E419FF", "#ECE51BFF", "#F5E61FFF","darkorchid1", "darkorchid2", "darkorchid3", "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF","#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#FF6699","#F68D45FF", "#FCA537FF","#F6E726FF", "#F4ED27FF", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#00489C", "#CCCCCC", "#999999","#999999", "#A1C299", "#300018","#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018")) + guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) + ylab("Relative Abundance (Family > 1%) \n") + xlab("Sample ID") + theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=8, colour="black"), axis.text.y = element_text(size=8, colour="black"), plot.title = element_text(hjust = 0.5), axis.ticks.x = element_line(colour="#000000", size=0.1), axis.ticks.y = element_line(colour="#000000", size=0.1), legend.text = element_text(size = 7))
print(PRSeqR_Season_family_bar)
p3 <- ggplotly(PRSeqR_Season_family_bar)
phtmlwidgets::saveWidget(p, "Season_family.html")
htmltools::tags$iframe(
src=file.path(getwd(), "Season_family.html"),
width="100%",
height="600",
scrolling="no",
seamless="seamless",
frameBorder="0")Creating Season vs. Location barplot with only normal samples
PRSeqR_family_Location <- PR_Normal %>%
tax_glom(taxrank = "Family") %>%
#Agglomerate at family level
transform_sample_counts(function(x) {x/sum(x)}) %>%
#Transform to rel.abundance
psmelt() %>%
# Melt to long format
group_by(Season, Location, Kingdom, Phylum, Class, Order, Family) %>%
dplyr::summarize(Mean = mean(Abundance, na.rm=TRUE)) %>% #Calculate average
filter(Mean > 0.01) %>%
#Filter arrange
arrange(Class) ## `summarise()` has grouped output by 'Season', 'Location', 'Kingdom', 'Phylum',
## 'Class', 'Order'. You can override using the `.groups` argument.#Creating a new column that combines both the Class and Family level information #so that Family level identifiers that are not unique to one Class (such as Ambiguous_taxa) aren't merged into a single large category when graphed
PRSeqR_family_Location$Class_family <- paste(PRSeqR_family_Location$Class, PRSeqR_family_Location$Family, sep="_") #Creating plot
PRSeqR_Location_family_bar <- ggplot(PRSeqR_family_Location, aes(x = Location, y = Mean, fill = Class_family)) + geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") + scale_fill_manual(values = c("black", "#440154FF", "#450659FF","#460B5EFF", "#472D7AFF", "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF", "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF", "#2E6E8EFF", "#2D718EFF", "#2B748EFF", "#2A778EFF", "#297B8EFF", "#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF", "#1FA287FF", "#21A585FF", "#23A983FF", "#25AC82FF", "#29AF7FFF", "#2DB27DFF", "#32B67AFF", "#37B878FF", "#3CBC74FF", "#57C766FF", "#5EC962FF", "#A2DA37FF", "#DAE319FF", "#E4E419FF", "#ECE51BFF", "#F5E61FFF", "darkorchid1", "darkorchid2", "darkorchid3", "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF", "#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#00CCCC", "#006600", "#FF6699", "#F68D45FF", "#FCA537FF", "#F6E726FF", "#F4ED27FF", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#191970","#353935","#F0FFFF","#89CFF0","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58")) + guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) + ylab("Relative Abundance (Family > 1%) \n") + xlab("Location") + theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=8, colour="black"), axis.text.y = element_text(size=8, colour="black"), plot.title = element_text(hjust = 0.5), axis.ticks.x = element_line(colour="#000000", size=0.1), axis.ticks.y = element_line(colour="#000000", size=0.1), legend.text = element_text(size = 7)) + facet_grid(. ~ Location, margins = TRUE, scale="free")
print(PRSeqR_Location_family_bar)
p4 <- ggplotly(PRSeqR_Location_family_bar)htmlwidgets::saveWidget(p,"PRSeqR_Location_family_bar")
htmltools::tags$iframe(
src=file.path(getwd(),"PRSeqR_Location_family_bar"),
width="100%",
height="600",
scrolling="no",
seamless="seamless",
frameBorder="0")Creating location vs. year barplot
PR_Normal <- PRphyseq %>% subset_samples(Sample_Type %in% c("N"))
PR_2022 <- PRphyseq %>% subset_samples(Year %in% c("2022"))
##Creating a subset of samples for only normal samples.This will exclude the enriched samplesPRSeqR_family_Year <- PRphyseq %>%
tax_glom(taxrank = "Family") %>%
transform_sample_counts(function(x) {x/sum(x)}) %>%
psmelt() %>%
group_by(Year,Season, Kingdom, Phylum, Class, Order, Family)%>%
dplyr::summarize(Mean =
mean(Abundance, na.rm=TRUE)) %>%
filter(Mean > 0.01) %>%
arrange(Class) ## `summarise()` has grouped output by 'Year', 'Season', 'Kingdom', 'Phylum',
## 'Class', 'Order'. You can override using the `.groups` argument.## `summarise()` has grouped output by 'Season', 'Kingdom', 'Phylum', 'Class',
## 'Order'. You can override using the `.groups` argument.#Creating a new column that combines both the Class and Family level information #so that Family level identifiers that are not unique to one Class (such as Ambiguous_taxa) aren't merged into a single large category when graphed
PRSeqR_family_Year$Class_family <- paste(PRSeqR_family_Year$Class, PRSeqR_family_Year$Family, sep="_") #Creating plot
PRSeqR_Year_family_bar <- ggplot(PRSeqR_family_Year, aes(x = Year, y = Mean, fill = Class_family)) + geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") + scale_fill_manual(values = c("black", "#440154FF", "#450659FF","#460B5EFF", "#472D7AFF", "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF", "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF", "#2E6E8EFF", "#2D718EFF", "#2B748EFF", "#2A778EFF", "#297B8EFF", "#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF", "#1FA287FF", "#21A585FF", "#23A983FF", "#25AC82FF", "#29AF7FFF", "#2DB27DFF", "#32B67AFF", "#37B878FF", "#3CBC74FF", "#57C766FF", "#5EC962FF", "#A2DA37FF", "#DAE319FF", "#E4E419FF", "#ECE51BFF", "#F5E61FFF", "darkorchid1", "darkorchid2", "darkorchid3", "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF", "#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#00CCCC", "#006600", "#FF6699", "#F68D45FF", "#FCA537FF", "#F6E726FF", "#F4ED27FF", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#191970","#353935","#F0FFFF","#89CFF0","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58")) + guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) + ylab("Relative Abundance (Family > 1%) \n") + xlab("Year") + theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=8, colour="black"), axis.text.y = element_text(size=8, colour="black"), plot.title = element_text(hjust = 0.5), axis.ticks.x = element_line(colour="#000000", size=0.1), axis.ticks.y = element_line(colour="#000000", size=0.1), legend.text = element_text(size = 7)) + facet_grid(. ~ Season, margins = TRUE, scale="free")
print(PRSeqR_Year_family_bar)
p5 <- ggplotly(PRSeqR_Year_family_bar)Note that 2023 is only enriched samples so this year cannot be accurately depicted and 2022 is missing from the wet season.
*Now that we know there is an error in the 2022 data when attempting to generate the 2022 NMDS plot, we will be creating a family plot for the 2022 samples to determine what is too similar - this is what we think is preventing the NMDS plot from being generated.
PRSeqR_2022_Sample <- PR_2022 %>%
tax_glom(taxrank = "Family") %>%
transform_sample_counts(function(x) {x/sum(x)}) %>%
psmelt() %>%
group_by(Location,Season,Kingdom,Phylum,Class,Order,Family,Genus,Species)%>%
dplyr::summarize(Mean =
mean(Abundance, na.rm=TRUE)) %>%
filter(Mean > 0.00) %>%
arrange(Class) ## `summarise()` has grouped output by 'Location', 'Season', 'Kingdom', 'Phylum',
## 'Class', 'Order', 'Family', 'Genus'. You can override using the `.groups`
## argument.##`summarise()` has grouped output by 'Season', 'Kingdom', 'Phylum', 'Class',
## 'Order'. You can override using the `.groups` argument#Creating a new column that combines both the Class and Family level information #so that Family level identifiers that are not unique to one Class (such as Ambiguous_taxa) aren't merged into a single large category when graphed
PRSeqR_2022_Sample$Family <- paste(PRSeqR_2022_Sample$Class, PRSeqR_2022_Sample$Family, sep="_") #Creating plot
PRSeqR_2022_Sample_bar <- ggplot(PRSeqR_2022_Sample, aes(x = Location, y = Mean, fill = Family)) + geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") + scale_fill_manual(values = c("black", "#440154FF", "#450659FF","#460B5EFF", "#472D7AFF", "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF", "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF", "#2E6E8EFF", "#2D718EFF", "#2B748EFF", "#2A778EFF", "#297B8EFF", "#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF", "#1FA287FF", "#21A585FF", "#23A983FF", "#25AC82FF", "#29AF7FFF", "#2DB27DFF", "#32B67AFF", "#37B878FF", "#3CBC74FF", "#57C766FF", "#5EC962FF", "#A2DA37FF", "#DAE319FF", "#E4E419FF", "#ECE51BFF", "#F5E61FFF", "darkorchid1", "darkorchid2", "darkorchid3", "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF", "#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#00CCCC", "#006600", "#FF6699", "#F68D45FF", "#FCA537FF", "#F6E726FF", "#F4ED27FF", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#191970","#353935","#F0FFFF","#89CFF0","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58")) + guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) + ylab("Relative Abundance (Family > 1%) \n") + xlab("Location") + theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=8, colour="black"), axis.text.y = element_text(size=8, colour="black"), plot.title = element_text(hjust = 0.5), axis.ticks.x = element_line(colour="#000000", size=0.1), axis.ticks.y = element_line(colour="#000000", size=0.1), legend.text = element_text(size = 7))
print(PRSeqR_2022_Sample_bar)
p6 <- ggplotly(PRSeqR_2022_Sample_bar)This is showing the samples averaged together, but there are multiple samples from each site. We are going to plot based on sample name to give us individual bars for each sample to better understand the similarities between each other.
PRSeqR_2022_Sample_ID <- PR_2022 %>%
tax_glom(taxrank = "Family") %>%
transform_sample_counts(function(x) {x/sum(x)}) %>%
psmelt() %>%
group_by(Sample_ID,Location,Season,Kingdom,Phylum,Class,Order,Family,Genus,Species)%>%
dplyr::summarize(Mean =
mean(Abundance, na.rm=TRUE)) %>%
filter(Mean > 0.00) %>%
arrange(Class) ## `summarise()` has grouped output by 'Sample_ID', 'Location', 'Season',
## 'Kingdom', 'Phylum', 'Class', 'Order', 'Family', 'Genus'. You can override
## using the `.groups` argument.## `summarise()` has grouped output by 'Season', 'Kingdom', 'Phylum', 'Class',
## 'Order'. You can override using the `.groups` argument#Creating a new column that combines both the Class and Family level information #so that Family level identifiers that are not unique to one Class (such as Ambiguous_taxa) aren't merged into a single large category when graphed
PRSeqR_2022_Sample_ID$Family <- paste(PRSeqR_2022_Sample_ID$Class, PRSeqR_2022_Sample_ID$Family, sep="_") #Creating plot
PRSeqR_2022_Sample_ID_bar <- ggplot(PRSeqR_2022_Sample_ID, aes(x = Sample_ID, y = Mean, fill = Family)) + geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") + scale_fill_manual(values = c("black", "#440154FF", "#450659FF","#460B5EFF", "#472D7AFF", "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF", "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF", "#2E6E8EFF", "#2D718EFF", "#2B748EFF", "#2A778EFF", "#297B8EFF", "#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF", "#1FA287FF", "#21A585FF", "#23A983FF", "#25AC82FF", "#29AF7FFF", "#2DB27DFF", "#32B67AFF", "#37B878FF", "#3CBC74FF", "#57C766FF", "#5EC962FF", "#A2DA37FF", "#DAE319FF", "#E4E419FF", "#ECE51BFF", "#F5E61FFF", "darkorchid1", "darkorchid2", "darkorchid3", "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF", "#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#00CCCC", "#006600", "#FF6699", "#F68D45FF", "#FCA537FF", "#F6E726FF", "#F4ED27FF", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#191970","#353935","#F0FFFF","#89CFF0","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58")) + guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) + ylab("Relative Abundance (Family > 1%) \n") + xlab("Sample_ID") + theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=8, colour="black"), axis.text.y = element_text(size=8, colour="black"), plot.title = element_text(hjust = 0.5), axis.ticks.x = element_line(colour="#000000", size=0.1), axis.ticks.y = element_line(colour="#000000", size=0.1), legend.text = element_text(size = 7))
print(PRSeqR_2022_Sample_ID_bar)
p7 <- ggplotly(PRSeqR_2022_Sample_ID_bar)htmlwidgets::saveWidget(p,"PRSeqR_2022_Sample_ID_bar")
htmltools::tags$iframe(
src=file.path(getwd(),"PRSeqR_2022_Sample_ID_bar"),
width="100%",
height="600",
scrolling="no",
seamless="seamless",
frameBorder="0")Creating barplots for enriched samples
#PR_ENR <- PRphyseq %>% subset_samples(Sample_Type %in% c("E"))
##Creating a subset of samples for only enriched samples.This will exclude the normal samples#PRSeqR_family_Location <- PR_ENR %>%
# tax_glom(taxrank = "Species") %>%
# transform_sample_counts(function(x) {x/sum(x)}) %>%
# psmelt() %>%
# group_by(Location,Season,Kingdom,Phylum,Class,Order,Family,Genus,Species)%>%
# dplyr::summarize(Mean =
# mean(Abundance, na.rm=TRUE)) %>%
# filter(Mean > 0.00) %>%
# arrange(Class)
## `summarise()` has grouped output by 'Season', 'Kingdom', 'Phylum', 'Class',
## 'Order'. You can override using the `.groups` argument#Creating a new column that combines both the Class and Family level information #so that Family level identifiers that are not unique to one Class (such as Ambiguous_taxa) aren't merged into a single large category when graphed
PRSeqR_family_Location$Genus_species <- paste(PRSeqR_family_Location$Class, PRSeqR_family_Location$Species, sep="_") #Creating plot ## Warning: Unknown or uninitialised column: `Species`.PRSeqR_Location_family_bar <- ggplot(PRSeqR_family_Location, aes(x = Location, y = Mean, fill = Genus_species)) + geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") + scale_fill_manual(values = c("black", "#440154FF", "#450659FF","#460B5EFF", "#472D7AFF", "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF", "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF", "#2E6E8EFF", "#2D718EFF", "#2B748EFF", "#2A778EFF", "#297B8EFF", "#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF", "#1FA287FF", "#21A585FF", "#23A983FF", "#25AC82FF", "#29AF7FFF", "#2DB27DFF", "#32B67AFF", "#37B878FF", "#3CBC74FF", "#57C766FF", "#5EC962FF", "#A2DA37FF", "#DAE319FF", "#E4E419FF", "#ECE51BFF", "#F5E61FFF", "darkorchid1", "darkorchid2", "darkorchid3", "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF", "#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#00CCCC", "#006600", "#FF6699", "#F68D45FF", "#FCA537FF", "#F6E726FF", "#F4ED27FF", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#191970","#353935","#F0FFFF","#89CFF0","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58")) + guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) + ylab("Relative Abundance (Species > 0) \n") + xlab("Location") + theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=8, colour="black"), axis.text.y = element_text(size=8, colour="black"), plot.title = element_text(hjust = 0.5), axis.ticks.x = element_line(colour="#000000", size=0.1), axis.ticks.y = element_line(colour="#000000", size=0.1), legend.text = element_text(size = 7))
print(PRSeqR_Location_family_bar)
p8 <- ggplotly(PRSeqR_Location_family_bar)#(2/2) This step is needed to create an interactive plot on gtihub. After it has been pushed to github you can download the html doc and it will open as the interactive plot via your web browser.
htmlwidgets::saveWidget(p, "PRSeqR_Location_family_bar")
htmltools::tags$iframe(
src=file.path(getwd(), "PRSeqR_Location_family_bar"),
width="100%",
height="600",
scrolling="no",
seamless="seamless",
frameBorder="0"
)Transforming Phyloseq object to relative abundance
#Transform to relative abundance
PR_RA_Normal = transform_sample_counts(PR_Normal, function(x){x / sum(x)})
#View the otu_table to see the transformed relative abundances
View(otu_table(PR_RA_Normal)) NMDS Plots
Choose colors and shapes
#Choose colors
plot.colors <- c("steelblue2","purple4","darkorange","firebrick","springgreen4", "gold", "darkblue", "yellowgreen","turquoise4", "orange","indianred","darkslategrey", "lightblue","darkgreen","mediumaquamarine","gray48","mediumorchid1", "#5F7FC7","#DA5724", "#508578", "#CBD588","#CD9BCD","#AD6F3B", "#673770","#D14285", "#652926", "#C84248", "#8569D5", "#5E738F","#D1A33D","#8A7C64", "#599861","dodgerblue","darkmagenta", "forestgreen","steelblue1", "cyan","mediumorchid3", "cadetblue3", "yellow", "#00FF7F","#40E0D0","#E97451","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58", "#C4B454")
plot.shape <- c(21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21)Subset the samples
Subsetting from PR_physeq_RA_Normal, which includes only the non-enriched samples. This dataframe has been filtered so that ASVs that were 0 in all samples are excluded.
#PR_S003 <- PR_physeq_RA_Normal %>% subset_samples(Project %in% c("S003"))
#PR_S015 <- PR_physeq_RA_Normal %>% subset_samples(Project %in% c("S015"))
#PR_2021 <- PR_Normal %>% subset_samples(Year %in% c("2021"))
#PR_2022 <- PR_Normal %>% subset_samples(Year %in% c("2022"))NMDS of both runs
all.nmds.source.ord <- ordinate(
physeq = PR_RA_Normal,
method = "NMDS",
distance = "bray")## Run 0 stress 0.1850681
## Run 1 stress 0.1850731
## ... Procrustes: rmse 0.0007928043 max resid 0.00578824
## ... Similar to previous best
## Run 2 stress 0.2139555
## Run 3 stress 0.2207744
## Run 4 stress 0.2022847
## Run 5 stress 0.2123585
## Run 6 stress 0.1850805
## ... Procrustes: rmse 0.001712935 max resid 0.01392611
## Run 7 stress 0.1871644
## Run 8 stress 0.1871644
## Run 9 stress 0.2024291
## Run 10 stress 0.2035301
## Run 11 stress 0.2080844
## Run 12 stress 0.1850833
## ... Procrustes: rmse 0.003884105 max resid 0.03295702
## Run 13 stress 0.1850062
## ... New best solution
## ... Procrustes: rmse 0.005042124 max resid 0.03300655
## Run 14 stress 0.2179661
## Run 15 stress 0.2297205
## Run 16 stress 0.1849911
## ... New best solution
## ... Procrustes: rmse 0.002372687 max resid 0.01513598
## Run 17 stress 0.1849877
## ... New best solution
## ... Procrustes: rmse 0.001707955 max resid 0.015439
## Run 18 stress 0.1850491
## ... Procrustes: rmse 0.005735748 max resid 0.05125967
## Run 19 stress 0.2237442
## Run 20 stress 0.1870003
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 1: no. of iterations >= maxit
## 18: stress ratio > sratmax
## 1: scale factor of the gradient < sfgrminNMDs plot of location for combined runs
#Plot, color coding by particle type
all.nmds.Location <- plot_ordination(
physeq = PR_RA_Normal,
ordination = all.nmds.source.ord) +
scale_fill_manual(values = plot.colors, "Location")+
scale_shape_manual(values = plot.shape, name = "Location") +
geom_point(mapping = aes(fill = factor(Location), shape = Location, size = 5)) +
guides(size=FALSE) +
guides(shape = guide_legend(override.aes = list(size = 3))) +
theme(plot.title = element_text(size = 18),
text = element_text(size = 18),
axis.title = element_text(size = 15),
panel.spacing = unit(1, "lines"),
panel.border = element_rect(colour = "black", fill = NA, size = 0.5),
panel.background = element_blank(),
legend.text = element_text(size = 15),
legend.title = element_text(size = 15),
legend.position = c(.99, .99),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
stat_ellipse(aes(group=Location))## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.print(all.nmds.Location)
NMDS plot of season
#Plot, color coding by seaosn
all.nmds.Season <- plot_ordination(
physeq = PR_RA_Normal,
ordination = all.nmds.source.ord) +
scale_fill_manual(values = plot.colors, "Season") +
scale_shape_manual(values = c("D" = 21, "W" = 21), name = "Season") +
geom_point(mapping = aes(fill = factor(Season), shape = Season, size = 5)) +
guides(size=FALSE) +
guides(shape = guide_legend(override.aes = list(size = 3))) +
theme(plot.title = element_text(size = 18),
text = element_text(size = 18),
axis.title = element_text(size = 15),
panel.spacing = unit(1, "lines"),
panel.border = element_rect(colour = "black", fill = NA, size = 0.5),
panel.background = element_blank(),
legend.text = element_text(size = 15),
legend.title = element_text(size = 15),
legend.position = c(.99, .99),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
stat_ellipse(aes(group=Season))
print(all.nmds.Season)
PERMANOVA Analysis
#filt_complete_bray <- phyloseq::distance(PR_RA_Normal, method = "bray")
#filt_complete_bray_df <- data.frame(sample_data(PR_RA_Normal))
#filt_complete <- adonis2(filt_complete_bray ~ mm_Hg*DO_mg_L*DO_LDO*Specific_Conductance_uS_cm, data = filt_complete_bray_df)
#filt_complete#Run PERMANOVA with Bray-Curtis dissimilarity
#filt_complete_bray <- phyloseq::distance(PR_RA_Normal, method = "bray")
#filt_complete_bray_df <- data.frame(sample_data(PR_RA_Normal))
#filt_complete <- adonis2(filt_complete_bray ~ Temperature_C*Location*pH*Salinity_ppt, data = filt_complete_bray_df)
#filt_complete#Run PERMANOVA with Bray-Curtis dissimilarity
#filt_complete_bray <- phyloseq::distance(PR_RA_Normal, method = "bray")
#filt_complete_bray_df <- data.frame(sample_data(PR_RA_Normal))
#filt_complete <- adonis2(filt_complete_bray ~ Temperature_C*Location*DO_LDO*pH*Salinity_ppt, data = filt_complete_bray_df)
#filt_complete#Convert output to dataframe
#filt.complete.pvalue <- as.data.frame(filt_complete)
#Create table of p-values
#knitr::kable(filt.complete.pvalue, row.names = NA) %>%
# kableExtra::kable_styling("striped",
# latex_options="scale_down") %>%
#kableExtra::scroll_box(width = "100%")PERMANOVA analysis of location vs pH
#Subset phyloseq object
#filt.mp.particle.RA <- subset_samples(MPfiltRA, sample_type == "Particle")
#sample_data(filt.mp.particle.RA)#Run PERMANOVA with Bray-Curtis dissimilarity
#filt_particle_bray <- phyloseq::distance(filt.mp.particle.RA, method = "bray")
#filt_particle_bray_df <- data.frame(sample_data(filt.mp.particle.RA))
#filt_all <- adonis2(filt_particle_bray ~ particle_type*effluent*week*polymer_type, data = filt_particle_bray_df)
#filt_all#Convert output to dataframe
#filt.all.pvalue <- as.data.frame(filt_all)
#Create table of p-values
#knitr::kable(filt.all.pvalue, row.names = NA) %>%
#kableExtra::kable_styling("striped",
# latex_options="scale_down") %>%
#kableExtra::scroll_box(width = "100%")#Run PERMANOVA with Bray-Curtis dissimilarity, only looking at particle type (glass vs. plastic)
#filt_particle_only_bray <- phyloseq::distance(filt.mp.particle.RA, method = "bray")
#filt_particle_only_bray_df <- data.frame(sample_data(filt.mp.particle.RA))
#filt_only_particle <- adonis2(filt_particle_only_bray ~ particle_type, data = filt_particle_only_bray_df)
#filt_only_particleNMDS Plot extras in R: Envfit
##Load libraries and read in your data
library(vegan)
library(ggplot2)##Subset your data to only environmental data and abundance data.
##I already had a dataframe from the PR_RA_Normal for abundance data, so I just needed to do the same for the sample data (aka environmental data).
#df_env <- as.data.frame(sample_data(PR_RA_Normal))
#df_com <- as.data.frame(otu_table(PR_RA_Normal))
## Individual data frames created, now I am selecting the columns that will be used for this plot and renaming the data frame to com and env.
#com = df_com[,1:95]
#env = df_env[,15:17, 20, 21,24]#convert com to a matrix and perform the NMDs ordination
#m_df_com = as.matrix(com)
#nmds code
#set.seed(123)
#nmds = metaMDS(m_df_com, distance = "bray")
#nmds##Running envfit function with env dataframe
#en = envfit(nmds,env, permutations = 999, na.rm = TRUE)RDA - Redundancy Analysis
library(phyloseq)
library(ggplot2)
library(microViz)
# install.packages("remotes")
#remotes::install_github("statdivlab/corncob")
#library(corncob)
#> microViz version 0.12.0 - Copyright (C) 2023 David Barnett
#> ! Website: https://david-barnett.github.io/microViz
#> ✔ Useful? For citation details, run: `citation("microViz")`
#> ✖ Silence? `suppressPackageStartupMessages(library(microViz))`
#knitr::opts_chunk$set(fig.width = 7, fig.height = 6)#PR_RA_Normal <- corncob::PR_phylo
#PR_phylo
#> phyloseq-class experiment-level object
#> otu_table() OTU Table: [ 36349 taxa and 91 samples ]
#> sample_data() Sample Data: [ 91 samples by 15 sample variables ]
#> tax_table() Taxonomy Table: [ 36349 taxa by 7 taxonomic ranks ]##Redundancy analysis is a constrained ordination method. It displays the microbial variation that can also be explained by selected constraint variables. Behind the scenes, a linear regression model is created for each microbial abundance variable (using the constraints as the explanatory variables) and a PCA is performed using the fitted values of the microbial abundances.
#PR_RA_Normal %>%
# ps_mutate(
# IBD = as.numeric(ibd == "ibd"),
# Female = as.numeric(gender == "female"),
# ) %>%
# tax_transform("clr", rank = "Genus") %>%
# ord_calc(
#constraints = c("IBD", "Female", "Abx."),
# method = "RDA", # Note: you can specify RDA explicitly, and it is good practice to do so, but microViz can guess automatically that you want an RDA here (helpful if you don't remember the name?)
# scale_cc = FALSE # doesn't make a difference
# ) %>%
# ord_plot(
# colour = "DiseaseState", size = 2, alpha = 0.5, shape = "active",
# plot_taxa = 1:8
# )